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Abstract — This letter considers the problem of recovering a 
positive stream of Diracs on a sphere from its projection onto 
the space of low-degree spherical harmonics, namely, from its 
low-resolution version. We suggest recovering the Diracs via a 
tractable convex optimization problem. The resulting recovery 
error is proportional to the noise level and depends on the density 
of the Diracs. We validate the theory by numerical experiments. 


I. Introduction 

Many applications in engineering and physics consider 
signals that lie on spheres (see for instance, a, m, 1271, 
l22l). In this letter we consider the problem of recovering a 
positive stream of Diracs on the sphere from its low-resolution 
measurement. The natural way to model a low-resolution 
version of a signal on a sphere is by its projection onto the 
space of low-degree spherical harmonics, as will be explained 
in Section In] 

This work is motivated by the problem of estimating the 
orientations of the white matter fibers in the brain using 
diffusion weighted magnetic resonance imaging (MRI) 1241 . 
It is common to model the measured signal as a spherical 
convolution of the underlying distribution of fiber bundles, 
called the orientation density function, with the point spread 
function of the diffusion tensor imaging sequence which 
smears out the fine details of the fibers’ distribution. The 
orientation density function is modeled as a stream of Diracs 
on the sphere. The locations and the positive weights of the 
Diracs represent the orientations of the fibers and the partial 
volume of the fiber within a voxel, respectively ED, ini. 
Therefore, the mathematical model elaborated in Section El 
suits this application. 

From the theoretical side, as far as we know, this is the first 
work to suggest a stable recovery of positive signals on the 
sphere from their low-resolution measurements. This result is 
part of an ongoing effort to derive recovery guarantees for 
super-resolution of signals in various geometries and settings 
(see, e.g. QOl, M, EQI, d, Ha, Q, 0, Q, 0). 

Several papers considered the recovery of Diracs on a sphere 
with general coefficients (not necessarily positive) from their 
low-resolution measurements. In 0, Q it was shown that 
recovery via convex optimization methods is robust under 
the assumption that the Diracs are sufficiently separated (see 

This work was funded by the European Union’s Horizon 2020 research 
and innovation programme under grant agreement ERC-BNYQ, by the Israel 
Science Foundation under Grant no. 335/14, and by ICore: the Israeli 
Excellence Center ’Circle of Light'. 


Theorem III.21) . The papers IfTSl . lfT9ll employ a finite rate of 
innovation framework to super-resolve Diracs on the sphere. 
This approach does not need any assumption on the Diracs’ 
distribution. However, these works have no robustness guar¬ 
antees. Additionally, in it was proven that separation is 
necessary for robust recovery in the presence of noise by any 
method. 

Following ll29l . we show that if the Diracs are known to 
be positive then the separation condition can be replaced by a 
weaker condition called Rayleigh regularity, which quantifies 
the density of the Diracs. We suggest recovering the Diracs 
via a tractable convex optimization problem. The resulting 
recovery error is proportional to the noise level and depends 
on the Rayleigh regularity of the signal. 

The letter is organized as follows. In Section|II]we formulate 
the problem and present necessary mathematical background. 
Section|ni]presents our main result, which is proved in Section 
m Section [V] shows some numerical experiments, which 
corroborate the theoretical results, and Section [VT] concludes 
the letter. 

II. Problem Formulation and Background 

Spherical harmonics play a key role in the analysis of sig¬ 
nals in a vast variety of tasks, analysis methods and sampling 
theorems, see for instance; IH, OTIl . llT4l . Il26l . Il25l . II 2 TII . Let 
Tn denote the space of homogeneous spherical harmonics 
of degree n, which is the resttiction to the bivariate unit sphere 
of the homogeneous harmonic polynomials of degree n in K^. 

Any point on the bivariate unit sphere is parametrized 
by X := {(j),0) G [0, 27r) x [0,7r]. The distance between two 
points Xi , Xj G is measured as 

p {xi, Xj) := arccos {xi • Xj). (ILl) 

Let 

Yji k — ^Pn^k (cOS 0) ^ k — 7Z, . . . , ri, 

be an orthonormal basis of Tn where Pn^k (x) is an 

associated Legendre polynomial of degree n and order k, and 


2n -I- 1 (n — |fc|)! 

47r (n + |fc|)! 

The functions {Y„ are the eigenfunctions of the Laplacian 
on S^, and thus can be understood as the extension of Fourier 
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analysis on the sphere. Any function g € L 2 (S^) can be 
expanded as ||3 


if X satisfies the separation condition of Definition \IL1\ then 
the solution f of 


gi^) — 'y ^ y ^ (ff) Yji^k) 

n—O k——n 


In this work, we consider a discrete positive signal of the 
form 


M 


/N = X! ’ Cm> 0, (n.2) 


where S [x] is the Kronecker delta function and X := {xm} is 
the signal’s support. We assume that the signal lies on some 
predefined grid C and that any pair of points on the 
grid Xi,Xj £ satisfy p{xi,Xj) >lfL for some L > I/tt. 
The higher L is, the larger the target resolution we want to 
achieve. 

The information we have on the signal is its projection onto 
the space of the low N spherical harmonics 

yn,k = ifjYn^k) + t]n,k, ( 113 ) 

n = 0,...,iV, k = —n,... ,n, 


where p := {gn,k} is some noise or model mismatch which 
is assumed to be bounded. In matrix notation we may write 


y = FNf + g^s:=F*j,y = Pm f + (HA) 

where y := Fm is a linear operator mapping a signal 

to its low N spherical harmonic coefficients, and the adjoint 
operator is given by F^z{x) = Y.n<N,\k\<n^ri,kYn,k{x). The 
operator Pm = F^Fm is the orthogonal projection onto the 
space of spherical harmonics of degree N, denoted by Vm- 
We aim to recover the sets {cm} , {xm} from the noisy low- 
resolution measurements (linT l. 

In recent papers Ih), Q, it was shown that signals of the 
form (III. 2b with general coefficients (namely, not necessarily 
positive values) can be recovered robustly from Vm by solving 
a tractable convex program. This holds provided that the 
signal’s support X satisfies the following separation condition; 

Definition II.l. A set of points A" C is said to satisfy the 
minimal separation condition if 


min p (xi-Xj) > 

Xi^Xj^X, i^j 


V 

N’ 


for a fixed separation constant v that does not depend on W, 
where p is defined in (III.lb . 


We will also make use of the notion of the super-resolution 
factor (SRF). The SRF quantifies the ratio between the reso¬ 
lution we want to achieve, specified by the grid spacing 1 /L, 
and the resolution we measure, namely, 

SRF := (II.5) 

The main result of a then states the following; 

Theorem II.2. Let X = {xm} C be the support of a 
signal of the form dmi with general coefficients Cm £ R. Let 
{yn,k} be as in ( 1/7. 3b with \\g\\(^^ < 5. For sufficiently large L, 


min || 5 |L subject to \\y - FmpWm < 

J- ^ 


satisfies 


/-/ 


< CoSRF^S, 


for some fixed constant Cq- 


( 11 . 6 ) 


Remark II.3. The minimal separation constant v was evaluated 
numerically to be 27r. The separation of ^ coincides with the 
spatial resolution of signals on spheres 1301 . 

The aim of this letter is to derive a stability result for the 
recovery of positive signals on the sphere from their low- 
resolution measurements. In this case, as presented in the next 
section, the separation condition can be replaced by a weaker 
condition called Rayleigh regularity. 


III. Main result 

In l|6l, it was proven that a positive signal on the sphere 
with M < N/2 (i.e. maximal cardinality of N/2) can be 
perfectly recovered from its noiseless projection onto Vm by 
solving a convex optimization problem. A general signal of 
cardinality M can be recovered from by an 

algebraic approach HD. However, both recovery results are 
not stable in the presence of noise. 

To derive a stability result for the positive case, we use 
the notion of Rayleigh regularity. A univariate signal with 
Rayleigh regularity r has at most r spikes within a resolution 
cell of size ^ for a separation constant p. In the multidi¬ 
mensional case, the definition of Rayleigh regularity is less 
intuitive (see discussion in IH) and can be interpreted as 
a density measure of the signal. For r = 1, the definition 
coincides with Definition III. II 

Definition III.l. We say that the set V C is 

Rayleigh-regular with parameters (/r, r ; Nffi) and write V £ 
N,L) if 

• V—Vi U • • • U 'Pr and Vi n Pj = 0 for all i j, 

• for all i = 1,..., r, Vi satisfies the separation condition 
of Definition HO with constant p. 

We denote the set of positive signals of the form (III.2b with 
support X £ TP'^^{p, r; N,L) as Tl+{p, r; N,L). 

The notion of Rayleigh regularity was used in ll29l to derive 
stability results for the recovery of positive signals from their 
low-degree Fourier coefficients. Our results can be seen as an 
extension to signals on spheres. 

In the sequel, we assume that the noise satisfies < 

S and suggest to recover the signal by solving the feasibility 
(convex) problem 

find 5 S subject to ||s — FjvffUfj < S, g >0. (III.l) 

Now, we are ready to present our main theorem. The theorem 
states that by solving the convex program (llll.lb . one can 
stably recover a positive signal on the sphere from its low- 
resolution measurements. The recovery error is proportional 
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to the noise level and depends on the Rayleigh regularity of 
the signal. 


Theorem III.2. Let f £ TZ+{i/r,r; N,L) be of the form 
( I//. 21 ) and consider the measurement model Then, for 

sufficiently large SRF, any solution f of MII.R satisfies 


f-f 




for some fixed constant Ci > 0. 


Corollary III.3. In the noiseless case, (5 = 0, the recovery is 
exact. 


As we show in the simulations, minimizing the £i norm in 
dnu) among all feasible solutions results in a low recovery 
error. 


IV. Proof of Theorem IIII. 2I 
The proof exploits the technique presented in ll29l (see also 
Bl). We commence by presenting the following Lemma which 
is a direct consequence of the construction in 151 : 

Lemma IV.l. Suppose that the set S = {^m} C satisfies 
the separation condition of Definition 1/7. 71 Then, for suffi¬ 
ciently large N there exists a polynomial q £ Vm obeying 
q{f) < 1 and 

Q (Cm) — 0, VCm & 

9(C) > CliVV(C,Cm)^ p(C,Cm)<a/7V, CmSS, 

9(C) > C 2 , i/'p(C,Cm) > cr/iV, VCm e S, 

for constants cr,Ci >0 and 0 < C 2 < L 

Set h := f — f C S| where / is the solution of the 
convex program (llll.lb . Let % ■.= £'B\ : h [x\ < O} and 

thus 77 C 27. By assumption / £ TZ+{vr,r] N,L) and con¬ 
sequently 77 £ 7Z^‘^^(ur,r; N,L). Therefore, by Definition 
HIO the set 77 can be presented as a disjoint union of r 
sets 77 = where 77» e W'^^{vr,l-N,L) for all 

i = 1,..., r. Observe that by simple rescaling N = N/r we 
have 

7^*‘^“(l/r, 1; N,L)= 1; N/r, L). 

Then, for each set 77i there exists an associated interpolating 
polynomial qi{x) £ V/v/r as given in Lemma BV. 1 1 

The key ingredient of the proof is the following construc¬ 
tion: 

r 

9 (a::) = n 

where a > 0 is a constant to be determined later. A product of 
spherical harmonics of degrees iVi, N 2 is a spherical harmonic 
of degree 7Vi -f N 2 and the computation of the corresponding 
representation is known as Clebsch-Gordan. Therefore, q £ 
Vn- We denote by qL\x\ the restriction of q{x) to the grid 
By construction, for all a: S 77 we have qi {x) = 0 for some 
i = 1,..., r. Therefore, 

r 

9l N = 9i (a) - Q! = -a. 

^We assume here that N/r is an integer for clarity. This assumption is not 
necessary for the results to hold. 


Additionally, for sufficiently large SRF (see ( III.5b ) we get that 
for all X £ 

9lN > -a 

= C/r-^’-SEF-^^ - a. 


By setting 

a := ^C/r-^^SEF-^^ < 1/2, (IV.l) 

we conclude that 


qL[x] = -a Vxe77, (IV.2) 

9l [x\ > a, Wx £ Si\77. 


Once we have constructed the appropriate polynomial q^, 
the rest of the proof follows directly. On the one hand, by 
dllBI l and (HIO i we get 


I(9l,MI = 


< 


\{PNqL, h)\ 
\{qL,PNh)\ 




Pn 



< 

< 


II9l||,^ 

26. 


PNf - S 



(IV.3) 


On the other hand, combining (IIV.2b with the fact that and 
h have the same sign pattern on we have 


\{QL,h)\ = 


^ qL[x]h[x] 


cGS? 


= (IV.4) 


> a\\h\\^^. 

Combining (IIV.3b . (IIV.4b and (IIV. lb we conclude that 

= iCf^r^^SEF^^S. 


V. Numerical experiments 

We now verify the theoretical results of this paper via 
numerical experiments. The convex optimization problems 
were solved using CVX ifTSll . In all experiments, we set the 
separation constant 1 / to he ^ and chose a uniform grid 

The signal support was generated as a union of r disjoint 
sets that were drawn randomly on the sphere, while keeping 
the separation requirements of Definition IIII.ll For each sup¬ 
port location, an associated amplitude was drawn randomly 
from a uniform distribution on the interval (0,10]. Then, we 
computed the projection of the signal onto Vn and added an 
iid normal additive noise. 

We solved the feasibility convex program ( llll.lb and chose 
the solution with minimal 7i norm over all feasible solutions. 
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(a) The underlying signal. (b) Noisy measurements. (c) The recovered signal. 

Figure V.l; An example for the recovery of a signal on the sphere from its projection onto Vn with parameters L = 60, N = 15, 
SRF=4, r = 3, M = 41 and SNR=30 db. 



4. 


VI. Conclusion 


20 30 40 50 60 70 60 90 100 


Figure V.2: The mean recovery error (over 10 experiments) 
as a function of the noise level with parameters N = 12, L = 
50 and r = 2. The blue asterisks and the red crosses represent 
the recovery error with and without minimizing the norm 
among all solutions of (IIII.ll i. respectively. 


In this letter, we proved that a discrete positive stream 
of Diracs on a sphere can be recovered robustly from its 
low-resolution measurements by solving a tractable convex 
optimization problem. The recovery error is proportional to 
the noise level and depends on the distribution of the Diracs 
on the sphere. 

In practice, signals that have sparse representation in a 
continuous dictionary might not have sparse representation 
after discretization ns. An obvious technique to alleviate this 
basis mismatch is by fine discretization, which will increase 
the recovery error significantly according to Theorem IIII.2I 
Proving a version of Theorem IIII.2I for continuous signals is 
therefore an important extension for future work. 


Rayleigh regulaiity 

r = 1 

r = 2 

r = Z 

r = 4 

Mean recovery error 

0.0026 

0.0148 

0.0285 

0.0584 

Max recovery error 

0.0059 

0.0365 

0.0452 

0.0699 


Table V.l; Mean recovery error (over 10 experiments for 
each value) of the solutions of (IIII.ll) with minimal norm 
as a function of the signal’s Rayleigh regularity with the 
parameters L = 50, IV = 12, SRF k 4 and SNR=30 db. 


Figure IV. II presents a recovery example of a signal with 
Rayleigh regularity of r = 3 and M = 41 in a noisy 
environment of SNR = iO db. In Figure |V2] we compare the 
output of the CVX program as a function of the noise level, 
with and without minimizing the £i norm among all feasible 
solutions. The recovery error was computed as the normalized 
error, i.e. 

1 ^ I . 

error =— ^ \f[q,p]- f[q,p] ■ 

p,q=l 

Table IV. II shows the mean recovery error as a function of the 
Rayleigh regularity parameter r with SNR = 30 and SRF = 
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